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Key words Abstract In many lichen-forming fungi, molecular phylogenetic analyses lead to the discovery of cryptic species 
chemotypes within traditional morphospecies. However, in some cases, molecular sequence data also questions the separation 
cryptic species of phenotypically characterised species. Here we apply an integrative taxonomy approach — including morphologi- 
haplotypes cal, chemical, molecular, and distributional characters — to re-assess species boundaries in a traditionally speciose 


group of hair lichens, Bryoria sect. Implexae. We sampled multilocus sequence and microsatellite data from 142 
specimens from a broad intercontinental distribution. Molecular data included DNA sequences of the standard fungal 
markers ITS, IGS, GAPDH, two newly tested loci (FRBi15 and FRBi16), and SSR frequencies from 18 microsatellite 
markers. Datasets were analysed with Bayesian and maximum likelihood phylogenetic reconstruction, phenogram 
reconstruction, STRUCTURE Bayesian clustering, principal coordinate analysis, haplotype network, and several dif- 
ferent species delimitation analyses (ABGD, PTP, GMYC, and DISSECT). Additionally, past population demography 
and divergence times are estimated. The different approaches to species recognition do not support the monophyly 
of the 11 currently accepted morphospecies, and rather suggest the reduction of these to four phylogenetic species. 
Moreover, three of these are relatively recent in origin and cryptic, including phenotypically and chemically variable 
specimens. Issues regarding the integration of an evolutionary perspective into taxonomic conclusions in species 
complexes, which have undergone recent diversification, are discussed. The four accepted species, all epitypified 
by sequenced material, are Bryoria fuscescens, B. glabra, B. kockiana, and B. pseudofuscescens. Ten species rank 
names are reduced to synonymy. In the absence of molecular data, they can be recorded as the B. fuscescens 
complex. Intraspecific phenotype plasticity and factors affecting the speciation of different morphospecies in this 
group of Bryoria are outlined. 
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INTRODUCTION of delimited clades (Jakob & Blattner 2006, Lumley & Sperling 


2011), or incomplete lineage sorting (Saag et al. 2014, Leavitt 


Accurate identification and characterization of species is the 
basis of communication, conservation, resources management, 
and material used in biological research. However, in groups 
of relatively recent origin, species delimitation is often difficult 
(Jakob & Blattner 2006, Leavitt et al. 2011, Lumley & Sperling 
2011). Organisms are always evolving, changing in response to 
either selective pressures or genetic drift, so that delimiting units 
to accord species names is not always clear (Naciri & Linder 
2015). Several phenomena can hinder species delimitation: 
phylogenetic/phenotypic mismatches (Articus et al. 2002, Mark 
et al. 2016, Pino-Bodas et al. 2016), ‘intermediate’ specimens 
between generally accepted taxa (Seymour et al. 2007), hybridi- 
zation (Konrad et al. 2002, Steinova et al. 2013), an absence 
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et al. 2016). Long-term reproductive isolation may produce 
structured, non-overlapping lineages, whereas an intraspecific 
phylogeny, as well as a recent or contemporary speciation 
event, may produce reticulated lineages (Abbott et al. 2016). 


The family Parmeliaceae is one of the most studied amongst 
lichenised fungi. It contains many genera with species delimita- 
tion problems, such as Cetraria aculeata (Lutsak et al. 2017); 
Letharia (Altermann et al. 2014), the Parmotrema reticulatum 
complex (Del-Prado et al. 2016), and Pseudephebe (Boluda et 
al. 2016). In some cases, a lack of correlation between geno- 
types and phenotypes has led to the recognition of cryptic 
species within morphologically indistinguishable or scarcely 
indistinguishable morphospecies (Molina et al. 2011a, b, Leavitt 
et al. 2012a, b, Singh et al. 2015, Boluda et al. 2016, Del-Prado 
et al. 2016), and so far, more than 80 cryptic lineages have been 
detected in Parmeliaceae (Crespo & Lumbsch 2010, Divakar et 
al. 2010). However, in other cases there is a mismatch between 
lineages revealed by standard DNA-barcoding markers and 
long-accepted morphospecies (Articus et al. 2002, Seymour 
et al. 2007, Velmala et al. 2014, Mark et al. 2016, Kirika et al. 
2016a, b, McMullin et al. 2016). 


In the morphologically similar ‘beard’ and ‘hair’ lichens of the 
Alectoria sarmentosa, Bryoria sect. Implexae, and Usnea 
barbata species complexes (Velmala et al. 2014, Mark et al. 
2016, McMullin et al. 2016), DNA sequences from standard 
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barcoding markers show that what were considered well delim- 
ited morphospecies are found admixed in a single lineage that 
may be interpreted as a single phylogenetic species. In such 
situations, many processes may be operative, including envi- 
ronmental plasticity (Boluda et al. 2016), hybridisation, ances- 
tral polymorphisms, incomplete lineage sorting (Joly et al. 
2009), limited value of neutral markers (Bekessy et al. 2003), 
or morphological variability mediated by low selective pressure, 
genetic drift, or huge population sizes (Hartl & Clark 2007). In 
these cases, the use of additional markers, especially highly 
variable ones (e.g., microsatellites), may contribute to an ex- 
planation of the underlying phenomena. 


Chemical characters, mainly the production of polyketides, were 
accorded major importance in species delimitation in lichen- 
forming fungi in the 1960s and 1970s (Hawksworth 1976, 
Lumbsch 1988). These compounds are formed by the fungal 
partner, and that expression can differ according to the position 
in a thallus or in pure culture. For almost 50 years, chemical 
products, generally linked to minor morphological differences, 
have been used to circumscribe species in Bryoria (Hawksworth 
1972, Brodo & Hawksworth 1977, Myllys et al. 2011, Velmala et 
al. 2014). The advent of molecular phylogenetics has enabled 
such species concepts to be tested, and they have proved par- 
ticularly wanting in one group of species, those placed in Bryoria 
sect. Implexae (Myllys et al. 2011, Velmala et al. 2014, Boluda 
et al. 2015). Velmala et al. (2014) provided DNA sequence data 
for 11 species in the section, and with the exception of B. glabra, 
all the other species were intermixed in clades with diverse, and 
not concordant, chemical and morphological features. Geneti- 
cally indistinguishable taxa (with the markers used), maintain 
distinctive phenotypes even when growing in physical contact 
with one another (Velmala et al. 2014, Boluda et al. 2015), so 
the variation cannot be attributed solely to ecological factors. 


A study on the morphospecies B. fuscescens in central Spain 
(Boluda et al. 2015) revealed specimens with the same nuclear 
internal transcribed spacer sequence (nulTS) but different ex- 
trolites (compounds formed on the surface of or excreted from 
hyphae). Subsequent fieldwork across Europe has revealed 
further combinations of extrolites, and also specimens sharing 
characters of additional morphospecies. In order to under- 
stand the evolutionary processes involved in B. fuscescens 
and related species we have adopted an integrative approach 
including morphological, distributional, and chemical data to- 
gether with DNA sequences from three standard loci (Schoch 
et al. 2012), two newly tested loci, and eighteen microsatellite 
(SSRs) markers (Nadyeina et al. 2014). We then analysed 
these datasets in a rigorous statistical framework to effectively 
integrate an evolutionary perspective into a revised and defen- 
sible taxonomic treatment. These studies are reported here, 
and we anticipate that the experience gained in this group of 
lichens will inform how other species complexes with similarly 
discordant datasets can be addressed. 


Table 2 Primer information used in Bryoria sect. Implexae. 


MATERIALS AND METHODS 


Sampling 

We examined 142 specimens from 14 countries in Europe, 
the Mediterranean Basin, and North and South America, rep- 
resenting 11 named morphospecies in Bryoria sect. Implexae 
(Table 1). Our dataset included 91 of the 97 specimens used 
by Velmala et al. (2014) in their revision of B. sect. Implexae. 
Newly obtained sequences are shown in bold in Table 1. Bryoria 
furcellata was used as outgroup to root the tree (Velmala et al. 
2014). Names used in the analyses follow the species concepts 
adopted in Velmala et al. (2014). 


Morphology and chemistry 


The newly studied specimens (Table 1) were examined mor- 
phologically under a Nikon SMZ-1000 dissecting microscope, 
and hand-cut sections studied with a Nikon Eclipse-80i com- 
pound microscope equipped with bright field and differential 
interference contrast (DIC). Habit photographs were taken 
with a Nikon 105 mm f/2.8D AF Micro-Nikkor Lens coupled to 
a Nikon D90 camera with daylight. Spot tests (K, C, and PD) 
and TLC were carried out following Orange et al. (2010). Sol- 
vent system C (200 ml toluene / 30 mL acetic acid) was used 
for TLC, with concentrated acetone extracts at 50 “C spotted 
onto silica gel 60 F254 aluminium sheets (Merck, Darmstadt, 
Germany). Spotted sheets were dried for 10 min in an acetic 
acid atmosphere to maximize resolution. Segments from the 
same lichen branch were used for both TLC and DNA extrac- 
tion to avoid the possible risk of taking samples from mixed 
collections. Morphological and thin layer chromatographic (TLC) 
analyses of the samples used in Velmala et al. (2014; Table 1) 
were taken from that study. 


DNA dataset 


The molecular dataset comprised DNA sequences and SSRs 
frequencies. DNA extraction was performed with the DNeasy 
Plant Mini Kit (Qiagen, Barcelona, Spain), following the manu- 
facturer's instructions. 


Eighteen fungal-specific microsatellites markers (Bi01, Bi02, 
Bi03, Bi04, Bi05, Bi06, BiO7, Bi08, BiO9, Bi10, Bi11, Bi12, Bi13, 
Bi14, Bi15, Bi16, Bi18 and Bi19) were amplified following Na- 
dyeina et al. (2014) using fluorescently labelled primers. Frag- 
ment lengths were determined on an ABI PRISM® 3130 Genetic 
Analyser (Life Technologies, Carlsbad, CA, USA). Genotyping 
was performed using GeneScan-500 LIZ as the internal size 
standard and GeneMapper v. 3.7 (Applied Biosystems, Foster 
City, CA, USA). 

For DNA sequencing, five loci were selected (Table 2), three 
commonly used as standard markers in fungi (ITS, IGS, and 
GAPDH), which were also used in Velmala et al. (2014), and 
two microsatellite flanking regions tested here for the first time 


Marker Description Primer forward (5'-3') Source Primer reverse (5'-3") Source 

ITS Internal transcribed spacers ITS1-F: Gardes & Bruns (1993) ITS4: White et al. (1990) 
of the nuclear rDNA including the CTTGGTCATTTAGAGGAAGTAA TCCTCCGCTTATTGATATGC 
5.88 region 

IGS Intergenic spacer of the IGS12b: Printzen & Ekman (2002) SSU72R: Gargas & Taylor (1992) 
nuclear rDNA AGTCTGTGGATTAGTGGCCG TTGCTTAAACTTAGACATG 

GAPDH Glyceraldehyde 3-phosphate Gpd1-LM: Myllys et al. (2002) Gpd2-LM: Myllys et al. (2002) 
dehydrogenase gene partial sequence ATTGGCCGCATCGTCTTCCGCAA CCACTCGTTGTCGTACCA 

FRBi15 Flanking region of Bryoria sect. FRBi15f: This paper FRBi15r: This paper 
Implexae microsatellite marker 15 GTCATAAGGGTATCAATCC TGAAAAGGTTTGGTGACTC 

FRBI16 Flanking region of Bryoria sect. FRBi16f: This paper FRBi16r: This paper 
Implexae microsatellite marker 16 CGAGGTTTCAGGAAAGGGAA AGGAAGTGATGTCGAGGT 
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(FRBi15 and FRBi16). Microsatellite flanking regions are vari- 
able non-coding DNA fragments that can contain phylogeneti- 
cal signal through a neutral molecular evolution (Zardoya et 
al. 1996, Chatrou et al. 2009). To explore this possibility, the 
flanking regions of the 18 microsatellite markers were checked 
upstream and downstream in the 454 pyrosequencing contigs 
used for microsatellite searching in Nadyeina et al. (2014). 
The variability of each region was assessed with the number 
of variable sites in contigs supported by 2—16 copies. From the 
36 regions (two for each of the 18 microsatellites), the most 
variable flanking regions were in Bi15 and Bi16, and specific 
primers were designed for those loci (Table 2). 


New DNA sequences (Table 1) were obtained using polymerase 
chain reactions (PCRs) as follows: a reaction mixture of 25 uL, 
containing 12 pL sterile water, 9 uL JumpStartTM REDTaq 
ReadyMix PCR Reaction Mix (Sigma-Aldrich, St Louis, MI, 
USA), 1.25 uL of each primer (forward and reverse) at 10 uM, 
anda 1.5 uL DNA template. Cycling conditions for ITS, GAPDH, 
FRBi15, and FRBi16 were 2 min at 94 “C; 35 cycles of 30s 
at 94 °C; 30 s at 56 °C; 2 min at 72 °C; and a final extension 
of 5 min at 72 °C. For IGS, the cycling process was: 2 min at 
94 °C; 15 cycles of 30 s at 94 °C, 30 s at 55 °C (decreasing 
1 °C each cycle down to 40 °C), 2 min at 72 °C, then 35 cy- 
cles of 30 s at 94 °C, 30 s at 55 °C; 90 s at 72 °C, and a final 
extension of 5 min at 72 °C. PCR products were checked and 
quantified on 1 96 agarose gel stained with ethidium bromide 
and cleaned using Exonuclease | and FastAP Thermosensitive 
Alkaline Phosphatase (Thermo Fisher Scientific, Waltham, MA, 
USA) according to the manufacturer's instructions. Sequencing 
was performed with labelling using BigDye Terminator v. 3.1 Kit 
(Applied Biosystems) as follows: 25 cycles of 20 s at 96 °C, 5s 
at 50 °C, and 2 min at 60 °C. PCR products were cleaned-up 
with the BigDye XTerminator Purification Kit (Applied Biosys- 
tems) according to the manufacturer's instructions. Sequences 
were obtained in an ABI PRISM 3130 Genetic Analyser (Life 
Technologies) and manually adjusted using DNA Workbench 
v. 6 (CLC bio, Aarhus, Denmark) and MEGA5 (Tamura et al. 
2011). Newly generated sequences were deposited in GenBank 
(Table 1 in bold). 


Clustering methodologies 


Phenetic analyses 


Two presence/absence (1/0) matrices were constructed, one 
for the extrolites detected by TLC, and another with morphology 
and geography data (Appendix 1). Morphological characters scored 
comprised those traditionally used to separate morphospecies 
in the group: 

1. pale/dark thallus colour; 

2. branching angles (acute/obtuse/mixed); 

3. soralia (absent/fissural/tuberculate/both); and 

4 pseudocyphellae (conspicuous/inconspicuous). 


For distributions, Old World vs New World was used. The 
R package cluster (Maechler et al. 2013) was used to obtain 
the dissimilarity matrix, and then the pvclust package (Suzuki & 
Shimodaira 2006) was run to obtain a phenogram (Zamora et 
al. 2013). Multiscale bootstrap resampling with 10000 bootstrap 
(bp) replicates was used to obtain approximately unbiased (au) 
p-values for branch supports. Groups were considered as sup- 
ported when bp values exceeded 70 or au values exceeded 95. 


Phylogenetic tree 


Alignments for each locus were performed using MAFFT v. 7 
(http://mafft.cbrc.jp/alignment/server/; Katoh & Standley 2013) 
with the G-INS-i alignment algorithm, a ‘1PAM/K = 2’ scoring 
matrix, with an offset value of 0.1, and the remaining para- 
meters set as default. Alignments were deposited in TreeBASE 


under accession nos TB2:S20007 (ITS, IGS, and GAPDH), 
TB2:S20005 (FRBi15), and TB2:S20004 (FRBi16). RDP v. 4 
(Martin et al. 2010) was used to detect potential recombination 
events, through the methods RDP (Martin & Rybicki 2000), 
GENECONV (Padidam et al. 1999), Chimaera (Posada & Cran- 
dall 2001), Maxchi (Maynard-Smith 1992), Bootscan (Gibbs et 
al. 2000, Martin et al. 2005), SiScan (Weiller 1998, Gibbs et 
al. 2000), PhyIPro (Weiller 1998), and 3Seq (Boni et al. 2007). 
Partitionfinder (Lanfear et al. 2012) was used to detect possible 
intra-locus substitution model variability, resulting in the splitting 
of the ITS region into ITS1, 5.88, and ITS2, and coding each 
codon position separately in GAPDH. Models of DNA sequence 
evolution for each locus partition were selected with jModeltest 
v. 2.0 (Darriba et al. 2012), using the Akaike information criterion 
(AIC, Akaike 1974). The best-fit model of evolution obtained 
was: ITS1 = TIM2, 5.88 = K80, ITS2 = TIM2ef + G, IGS = 
TrN + 1, GAPDH jet position = TrN + I, GAPDH 2nd position = 
F81 + I, GAPDH 3th position = TPM3uf, FRBi15 = TPM3uf + I, 
FRBi16 = TPM3uf + G. To detect possible topological conflicts 
among loci, the CADM test (Legendre & Lapointe 2004, Camp- 
bell et al. 2011) was performed using the function ‘CADM.global’ 
implemented in the library ‘ape’ of R (Paradis et al. 2004). As 
loci FRBi15 and FRBi16 were not congruent among them and 
neither with the remaining loci, three alignments were used, 
resulting in three trees, one for each FRBi region and another for 
the concatenated dataset including loci ITS, IGS and GAPDH. 
For the concatenated matrix, specimens with more than one 
missing locus were excluded. Datasets were analysed using 
maximum likelihood (ML) and Bayesian (B/MCMCMC) ap- 
proaches with gaps treated as missing data. 


For ML tree reconstruction, we used RAXML v. 8.2.10 (Stama- 
takis 2006) implemented in CIPRES Science Gateway (https:// 
www.phylo.org/; Miller et al. 2010) with the GTRGAMMA 
model (Stamatakis 2006, 2014, Stamatakis et al. 2008). Sup- 
port values were assessed using the ‘rapid bootstrapping’ 
option with 1000 replicates. For the Bayesian reconstruc- 
tion, MrBayes v. 3.2.1 (Ronquist & Huelsenbeck 2003) was 
used. Two simultaneous runs with 10 M generations each, 
starting with a random tree and employing 12 simultaneous 
chains, were executed. Every 500th tree was saved to a file. 
Preliminary analysis resulted in an overestimation of branch 
lengths and to correct this we used the uniform compound 
Dirichlet prior brlenspr = unconstrained : gammadir (1, 1, 1, 1; 
Zamora et al. 2015). We plotted the log-likelihood scores of 
sample points against generations using Tracer v. 1.5 (Rambaut 
et al. 2014) and determined that stationarity had been achieved 
when the log-likelihood values of the sample points reached 
an equilibrium and ESS values exceeded 200 (Huelsenbeck 
& Ronquist 2001). Posterior probabilities (PPs) were obtained 
from the 50 96 majority rule consensus of sampled trees after 
excluding the initial 25 Yo as burn-in. The phylogenetic tree was 
drawn with FigTree v. 1.4 (Rambaut 2009). 


STRUCTURE 


STRUCTURE v. 2.3.4 (Pritchard et al. 2000, Falush et al. 2003) 
was run with the SSRs data matrix. Analysis was computed 
with 100000 burn-in generations and 100000 iterations using 
a K value from 1 to 12 (i.e., the putative number of species 
we may have) and 20 replicates for each K. To combine the 
20 runs of each K in a single result, CLUMMP v. 1.1.2 (Jakobs- 
son & Rosenberg 2007) was used and visualised replacing the 
CLUMMP output values in a STRUCTURE output of the same 
K, and then plotted using the STRUCTURE software. To show 
the probability of each K value, STRUCTURE HARVESTER 
(Earl & Von Holdt 2012), with the AK method (Evanno et al. 
2005) was used, considering the most probable K the first one 
that appears close to 0 in the output graphic. 
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Principal coordinate analysis 


Principal coordinate analysis (PCoA) was carried out with the 
SSRs length data in GenAIEx 6.5. The results of the three first 
axes were plotted in a three-axis graph using The Excel 3D 
Scatter Plot v. 2.1, in which the graphic can be moved in 3D to 
obtain a better understanding of how the plots are distributed in 
the space. Since the projection of this 3D graph on a paper is 
necessarily confusing, PCoA results were plotted on two differ- 
ent 2D graphs showing axes 1 and 2, and 1 and 3, respectively. 


Haplotype network 


Haplotype network reconstruction was performed using TCS 
v. 1.2.1 (Clement et al. 2000) with the concatenated sequences 
matrix, excluding the outgroup, using gaps as missing data, and 
a 95 % connection limit. Specimens differing only by missing or 
ambiguous characters were not counted as haplotypes. 


Species delimitation analyses 


In order to examine species delimitation, four computational 
approaches not requiring prior hypothesis of a putative number 
of species were used: 

1. Automatic Barcode Gap Discovery ABGD (Puillandre et al. 
2011) based on barcode gaps using genetic distances; 

2. Poisson Tree Processes PTP (Zhang et al. 2013), based 
on gene trees; 

3. The Generalized Mixed Yule coalescent approach GMYC, 
which combines a coalescent model of intraspecific branch- 
ing with a Yule model for interspecific branching (Pons et 
al. 2006, Monaghan et al. 2009); and 

4 DISSECT (Jones et al. 2014) based on the multispecies 
coalescent model for species delimitation. 


ABGD and PTP were carried out using the online servers http: A 
wwwabi.snv.jussieu.fr/public/abgd/ and http://species.h-its. 
org, respectively. GMYC was analysed with the gmyc func- 
tion in the SPLITS package in R (v. 2.10, www.cran.r-project. 
org), employing the single (GMYCs) and multiple (GMYCm) 
threshold methods. Because GMYC needs a strictly ultramet- 
ric and bifurcating tree with no zero branch lengths, identical 
sequences were deleted and an ultrametric tree was generated 
using BEAST v. 1.8.2 software (Drumond et al. 2012), with the 
evolutionary models explained in the Bayesian phylogenetic 
reconstruction. Arun of 100 M iterations logging every 1 000th 
iteration was conducted. Consensus tree was generated with 
TreeAnotator v. 1.8.2 after discarding the initial 10 96 trees as 
burn-in. ESS values above 200 were ensured using Tracer 
v. 1.6 (Rambaut et al. 2014). 


DISSECT analysis was implemented in STARBEAST (*BEAST, 
Drumond et al. 2012) using the concatenated DNA matrix after 
removing identical sequences and following the instructions of 
Jones et al. (2014). First, we used BEAUti (Drumond et al. 2012) 
to produce the xml file, with every individual encoded as if it was 
a separate species. Sites, clocks and trees were released as 
unlinked. Nucleotide substitution models and other parameters 
(as in the Bayesian analysis, see above), were encoded using 
BEAUti if possible, or manually entered. For the ITS locus, a 
substitution rate of 0.0033 substitutions per site per million years 
was introduced (Leavitt et al. 2012a, b), setting other loci as 
estimated with a lognormal relaxed clock. A birth-death-collapse 
prior that controlled the minimal split heights for the putative 
resulting species was manually added to the xml file. This prior 
contained the parameters CollapseHeight (£) with a value of 
0.0001 and CollapseWeight (w), set as estimated using a Beta 
distribution with values 10 and 1.5. Selected parameters provide 
a highest probability density around 4—5 clusters, the most 
probable number of taxa meriting separation according to other 
analyses performed for this paper. However, this prior is diffuse 
and allows to obtain a different number of putative taxa if they 


adjust better to the data. The xml file was executed in BEAST 
with 250 M MCMC iterations, sampling every 10 000th iteration. 
Tracer v. 1.6 (Rambaut et al. 2014) was used to assess ESS 
values above 140. The resulting *BEAST species tree output 
was then treated with SpeciesDelimitationAnalyzer (Jones et 
al. 2014), with a burn-in of 5000 trees (20 96 of the total gene- 
rated), a collapse height of 0.001 (one fraction lower than in 
the *BEAST analysis) and a simcutoff value of 1 to ignore this 
parameter, as according to sequence variability, we expected 
very similar putative species to emerge. The resulting similarity 
matrix was plotted with R v. 2.15.1 (R Core Team 2014) follow- 
ing the method of Jones et al. (2014). 


Divergence time estimation 


Two divergence time estimations were performed, one only with 
the ITS region and a defined substitution rate, and the other 
with the concatenated data matrix of ITS, IGS, and GAPDH 
loci. A rate of 3.30 x 10—9 s:s-1-yr—1 for the ITS region as a 
whole was used, with a GTR + G + | substitution model (Leavitt 
et al. 2012a, b). In the concatenated matrix analysis, as no 
previous literature on substitution rates for IGS and GAPDH in 
lichen-forming fungi is available, these were set as estimated 
in the Bayesian phylogenetic analysis. A *BEAST analysis was 
executed, using a relaxed clock model (uncorrelated lognormal), 
a birth-death model prior for the node heights and unlinked 
substitution models, clocks and trees for each partition. Clades 
G, Ko, NA, and WD were selected as potential species, forc- 
ing them to remain monophyletic (Fig. 6). No calibration points 
could be used, as no fossils or previous dating of this species 
complex are available. To avoid stochastic events, two inde- 
pendent analyses were run, each with 200 million generations, 
sampling each 5000 trees, and discarding the first 10000 trees 
(25 96) as burn-in. Tracer v. 1.6 (Rambaut et al. 2014) was used 
to ensure ESS parameter values above 115 in the concatenated 
matrix, and 185 for the ITS analysis. Different priors were 
tested but no higher ESS values could be obtained, which we 
suspect was due to the very similar sequences, and the uncer- 
tain topology of the backbone connecting the groups Ko, NA, 
and WD. The two runs performed for each input were merged 
with logcombiner v. 1.8.2, and the resulting trees merged in a 
consensus tree using TreeAnnotator v. 1.8.2 (Drumond et al. 
2012). FigTree v. 1.4 (Rambaut 2009) was used to display the 
ITS and the concatenated dated species trees. 


Demography 

Changes in population sizes through time were estimated us- 
ing the Bayesian skyline analysis (Drumond et al. 2005) with 
BEAST. Only clades Ko, NA, and WD, isolated and merged, 
were studied, as they show a clock-like tree topology and ad- 
equate sampling sizes. 


Following the methods used for divergence time estimation 
analysis, the demography analyses were run using the ITS 
region without partitioning, with the GTR + G + | model of 
nucleotide substitution and a substitution rate of 3.30 x 10—9 
s:s—1:yr—1, and with a strict molecular clock model (Leavitt et 
al. 2012a, b). Additionally, the same analysis was repeated with 
the concatenated data matrix using the ITS substitution rate, 
estimating the other loci rates with a relaxed clock model and 
using the nucleotide substitution models for IGS and GAPDH 
explained in the Bayesian phylogenetic reconstruction. Four 
independent runs for each input were processed with 50 M 
MCMC generations, sampling parameter values every 5000th 
generation, using the Bayesian Skyline tree prior model, six 
discreet changes in population size and the linear growth op- 
tion. ESS values were checked with Tracer v. 1.5 (Rambaut et 
al. 2014), and the two best of the four runs were combined, 
obtaining values usually above 200, with some exceptions with 
a lower limit of 100. Skyline plots were drawn with Tracer v. 1.5. 
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To support the Bayesian skyline test, a neutrality test was per- 
formed to infer if populations are in mutation-drift eguilibrium. 
Tajima's D (Tajima 1989) and Fu's Fs (Fu 1997) were calculated 
with DnaSP v. 5.10 (Librado & Rozas 2009). A significantly 
positive D is interpreted as a diversifying selection or a recent 
bottleneck, whereas a negative significant D shows purifying 
selection or a recent expansion. If D is not significantly different 
from 0, a mutation-drift equilibrium may be occurring. Fu's Fs 
can be interpreted in the same way. 


RESULTS 


Morphological and chemical clustering 


The wide geographical range of collections revealed a combi- 
nation of characters not previously reported in the Bryoria fus- 
cescens complex, especially those from the previously less- 
studied Mediterranean Basin. Specimens with intermediate 
morphologies amongst traditionally accepted species were rec- 
ognized, and the application of species names according to the 
current taxonomy was ambiguous. Individuals connecting the 
phenotypes and chemotypes of the taxa currently recognized as 
Bryoria fuscescens, B. implexa, B. kuemmerleana, and B. vran- 
giana were not rare in the Mediterranean Basin. For example, 
chemotypes thought to be diagnostic for a particular taxon 
were detected in specimens morphologically belonging to other 


FRBi15 


03 
|B. vrangiana S10 
IB. vrangiana LO5.17, 


ESBS BSED Po POP By EAD HOES ED TD SENSE ES 
3 3333 * 
H F SS E Š 1 
38055 SSeS EI 
EI ra Ei SS = 
EE SE 
E e: "Sep i 
à a 
S d DN 
2 Dx 
è 28s 
Soop OODD Bs 
nxn SLL) 2 
Bees D 
ZE 
SS 


. pikei pil 
B. vrangiana L03.07 
Bryona vrangiana L10.13 


0006 


taxa, as well as specimens producing extrolites characteristic 
of different taxa in a single thallus. Thin-layer chromatography 
(TLC) revealed seven different extrolites: alectorialic, barbatolic, 
fumarprotocetraric, gyrophoric, norstictic, and psoromic acids, 
atranorin, and sometimes also related substances such as 
chloroatranorin, protocetraric, or connorstictic acids. Atranorin, 
a typical accessory substance in the genus Bryoria, was not 
used in the posterior analyses because it appears in trace 
amounts in many samples and is often difficult to unequivocally 
discern if it is present or absent by TLC alone. 


The chemical presence/absence matrix resulted in the pheno- 
gram shown in Appendix 2a. The matrix included specimens 
with as many as four extrolites, something not previously 
reported in the complex. Chemical characters were separated 
into two main groups: 

1. specimens that contain benzyldepsides (i.e., alectorialic 
and barbatolic acids), substances traditionally used to se- 
parate B. capillaris and B. pikei from other species in the 
complex; and 

2. specimens without benzyldepsides. 


The latter were clustered in two well-supported groups, one with 
fumarprotocetraric acid as the main substance, and the other 
without it (including specimens with no detectable substances). 
If the structural relationships of the compounds were encoded in 
the presence/absence matrix (benzyldepsides vs depsidones), 
the same clustering was obtained. 
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Fig. 1 Phylogenetic relationships in Bryoria sect. Implexae based on FRBi15 and FRBi16 loci. Tree topology depicts the result of the Bayesian Markov chain 
Monte Carlo (B/MCMC) analysis. Posterior probabilities and bootstrap analysis for the supported nodes (2 0.95 and 2 70 %) are indicated at the main nodes. 
Lines connecting clades indicate putative recombination events, with main parents (continuous lines) and minor parents (discontinuous lines). Because the 
clade insertion in the trees is influenced by the recombination, clades with recombination are depicted with a discontinuous branch line. Note that clades with 
recombination appear as sister or close to the main parent but tending to be deviated towards the minor parent. — The coloured bar corresponds to the SSR 
genepool from Fig. 2, with specimens intermediate between two or more genepools in grey. The clades obtained, although well-supported, do not follow any 


evident geographic, morphologic, chemical, or microsatellite pattern. 
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The analysis of combined morphological, geographical, and 
chemical characters resulted in the phenogram in Appendix 2b. 
Only terminal branches were supported, including few mono- 
phyletic morphospecies, although not well isolated from others. 
Neither accepted morphospecies nor an unequivocal number 
of phenotypic groups could be recognized. This ambiguity was 
largely attributable to phenotypically intermediate specimens, 
mainly from the Mediterranean Basin, and also by the presence 
of some shared characters amongst the morphospecies, such 
as the presence/absence of soralia, pseudocyphellae, extrolite 
composition, and thallus colour. 


Phylogenetic tree 


Due to the topological conflict between loci, three DNA matrices 
were used to generate three phylogenetic trees: 

1. a concatenated matrix including ITS, IGS, and GAPDH 
with 134 individuals consisting of 1774 unambiguously 
aligned nucleotide position characters, with 83 parsimony 
informative (Pi) sites; 


2. FRBi15 with 93 individuals contained 569 unambiguously 
aligned nucleotide position characters, with 44 Pi sites; and 

3. FRBi16 with 80 individuals had 632 unambiguously aligned 
nucleotide position characters, with 160 Pi sites. 


No evidence of recombination events was detected in the 
concatenated matrix. The resulting tree (Fig. 6) had four well- 
supported main clades, G (Glabra, yellow), Ko (Kockiana, 
magenta), NA (North American, blue), and WD (Widely Dis- 
tributed, red + green + brown). Clade G included only material 
of B. glabra, appearing as an isolated taxon sister to the other 
three clades, which showed an uncertain topology between 
them. Clade Ko included material named B. kockiana and two 
unidentified specimens, all collected in Alaska (USA). Clade NA 
comprised the previously recognized North American morphos- 
pecies group (the ‘North American endemic species’, Velmala 
et al. 2014) named as B. friabilis, B. inactiva, B. pikei, and B. 
pseudofuscescens. While these species were mixed in the tree, 
the group as a whole was resolved as monophyletic. The WD 
clade included specimens widely distributed but mainly from 
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Fig.2 Bayesian inference of population structuring using STRUCTURE v. 2.3.4 (Pritchard et al. 2000, Falush et al. 2003) and nine microsatellite loci in Bryoria 
sect. Implexae. — Left. Results from the hypothesis of 2—6 clusters. Vertical bars represent specimen assignment probability into a genetic cluster depicted by 
the colours. Morphospecies names given to the specimens appear at the top. — Right. Detailed columns of the K - 6 hypothesis, the numbers representing the 
specimens shown in Table 1 to provide a better understanding of the components of each individual. — G = Glabra clade; Ko = Kockiana clade; NA = North 


American clade; WD = Widely Distributed clade; * = Bryoria pikei specimen 49. 
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Europe (the ‘European and globally distributed species’ group, 
Velmala et al. 2014) under the names B. capillaris, B. fusces- 
cens (syn. B. chalybeiformis and B. lanestris), B. implexa, B. 
kuemmerleana, and B. vrangiana. None of these previously 
recognised species formed a monophyletic group. 


The phylograms produced using the FRBi15 and FRBi16 mark- 
ers had a different tree topology, not congruent among them 
or with that from the preceding concatenated dataset. In the 
FRBi15 reconstruction (Fig. 1), B. glabra was not represented 
due to the lack of primer annealing in the PCR process, and 
the tree could not be rooted. Several well-supported groups 
were produced, but did not follow any evident geographic, mor- 
phological, chemical, or SSR frequency pattern. Bryoria pikei 
L376 had a sequence with a putative recombination fragment 
with the B. vrangiana S10 clade in c. 50 % of the total length. 
This insertion placed the specimen out of the main parental 
group and it appears as its sister. Although marker FRBi16 (Fig. 
1) also produced a well-resolved tree with supported nodes, 
the clades do not show any phenotypic and/or geographic 
structure. In this tree reconstruction, B. glabra was an isolated 
taxon and served to root the tree. In FRBi16 sequences, many 
putative recombination events were detected, suggesting a 
reticulate evolution. In both trees in Fig. 1, clade Ko (magenta) 
was recovered as monophyletic, but embedded between other 
named morphospecies. 


STRUCTURE clustering 


Of the 18 microsatellite markers, the nine that showed more 
than 95 Yo successful amplification across the samples were 
used (number of haplotypes shown in brackets, Appendix 3): 
Bi01 (17), BiO3 (6), Bi04 (8), BiO5 (5), Bi10 (8), Bi11 (9), Bi12 
(12), Bi14 (6), and Bi19 (5). We allowed a maximum of three 
missing loci per specimen, a value reached only in seven 
samples. STRUCTURE was allowed to run to K = 12, but from 
K = 6 the clustering process started to be uninformative (Fig. 2). 
The likelihood results of the AK analysis (Evanno et al. 2005) 
indicated three as the most probable number of clusters (likeli- 
hood = —1232, AK = 2.2), the clades G, NA, and WD (Fig. 2, 
K = 3). Clade Ko, which appeared isolated in the concatenated 
phylogenetic tree (Fig. 6), could not be accepted as distinct 
under a K = 3 hypothesis. However, B. glabra, a morphologi- 
cally delimited taxon, was not isolated at firstin STRUCTURE. 
This could be attributable to the clustering algorithms being 
influenced by unbalanced sampling sizes, masking clade Ko, 
which appeared isolated at K = 6. From K 2 4 to K= 6, the 
new groups appeared mainly inside the WD clade, showing 
that the samples from Europe were much more diverse than 
those from North America. Indeed, the NA clade was not split 
into subgroups even at K = 10. Apart from B. glabra, no other 
named morphospecies formed an exclusive cluster even at 
high K values. 
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Fig. 3 Principal Coordinate Analysis (PCoA) of microsatellite data in Bryoria sect. Implexae. Species names according to Velmala et al. 2014 (shape and 
colours) and STRUCTURE clusters (colours) for each specimen are represented in the three main coordinates. Note that the Ko clade does not appear iso- 
lated from the NA clade in any coordinate axis. — G = Glabra clade; Ko = Kockiana clade; NA = North American clade; WDr = Widely Distributed red clade; 


WDg = Widely Distributed green clade. 


Table3 Species delimitation analysis results for loci ITS, IGS, GAPDH and the concatenated data matrix in Bryoria sect. Implexae. Brackets indicate groups 
predicted as conspecific. — G = Glabra clade; Ko = Kockiana clade; NA = North American clade; WD = Wide Distributed clade; WDr = Wide Distributed red 


clade; WDg = Wide Distributed green clade; pik5 = Specimen Bryoria pikei 5. 


Method ITS IGS GAPDH Concatenated 

ABGD 2 =G + (Ko, NA, WD) 2 =G + (Ko, NA, WD) 4=G+Ko+NA+WD 4=G+Ko+NA+WD 

PTP 2 =G + (Ko, NA, WD) 2 =G + (Ko, NA, WD) 4=G+Ko+NA+WD 5 = G+ Ko + NA + WDr + WDg 
GMYCs 4 =G + (Ko, NA, WDg) + WDr + WDr 3 = G + (Ko, WD) + NA 4=G+Ko+NA+WD 6 =G + Ko + NA + pik5 + WDr + WDg 
GMYCm 4 =G + (Ko, NA, WDg) + WDr + WDr 4 = G + (Ko, WD) + NA1 + NA2 4=G+Ko+NA+WD 5 = G+ Ko + NA+ WDr + WDg 
DISSECT - - = 5= G + Ko + NA+ pik5 + WD 
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Principal coordinates analysis (PCoA) 


The PCoA analysis has a three-dimensional output represented 
here in two graphs, one comparing axis 1 against 2, and the 
other 1 against 3 (Fig. 3). The information percentage of each 
axis was 44.47 %, 15.06 %, and 14.44 %, respectively. Clade 
G (Fig. 3, yellow) appeared isolated, whereas clade Ko (Fig. 
3, magenta) was admixed with NA clade (Fig. 3, blue), form- 
ing a single cluster. Clade WD was isolated from the others, 
but divided into two clusters, one corresponding to the red 
and brown groups in the K = 6 STRUCTURE output (Fig. 2), 
and one for the green group. Apart from B. glabra, none of the 
currently accepted morphospecies formed a defined group. 
Four reasonably isolated clusters could be distinguished, cor- 
responding to the groups G, WDr (Widely Distributed, Fig. 3, 
red), WDg (Widely Distributed, Fig. 3, green), and Ko together 
with NA forming a single cluster. 


Haplotype network 


The haplotype network of the concatenated data matrix, cod- 
ing gaps as missing data, produced 39 haplotypes. Bryoria 
glabra specimens (Appendix 4, yellow) formed two haplotypes 
not connected with other members of the network, indicat- 
ing genetic isolation of this species. One of the haplotypes 
was composed exclusively of South American specimens, 
whereas the other contained European, North American, and 
South American samples. Clade Ko (Appendix 4, magenta) fell 
into two haplotypes, one including specimens with psoromic 
acid and identified as B. kockiana, and the other clustering 
unidentified samples with no substances detected. This group 
was connected to the NA clade (Appendix 4, blue) by a long 
branch with 13 mutation steps. The NA clade was separated 
by nine mutations from the WD clade (Appendix 4, green, red, 
and brown). The WD green, red, and brown groups split by 
STRUCTURE (Fig. 2) formed a unique cluster. Four isolated 
clusters could be distinguished, corresponding to the groups 
G, Ko, NA, and WD. 


Species delimitation programs 

The ABGD, PTP, GMYC, and DISSECT programs (Table 3) 
use different algorithms, and consequently different numbers 
of putative species may be predicted. The genetic distance 
method (ABGD) gave the smallest number of putative species, 


Fig. A Similarity matrix from DISSECT analysis performed after clone cor- 
rection in Bryoria sect. Implexae. Squares represent posterior probability 
(white = 0, black = 1) of pairs of specimens to belong to the same species. 
Resulting major groups are delimited by lines, which indicate the clade on 
the collapsed phylogenetic tree. 


whereas the coalescence methods (especially GMYC) the 
largest. Analyses also revealed the contribution of each locus 
to the postulated species delimitation, GAPDH being the most 
informative and constant marker. DISSECT analysis (Fig. 4) 
predicted five species corresponding to G, Ko, NA, and WD 
clades, and specimen B. pikei 5. Although the GMYC analysis 
also showed the B. pikei 5 specimen as a separate species, it 
was grouped in the NA clade in the other analyses. DISSECT 
showed two internal greyish square groups in WD, but they did 
not correspond exactly to the WDr and WDg groups in Fig. 2, 3 
(STRUCTURE and PCoA analyses). 


Node dates and demographic history 


The calibrated maximum clade credibility chronogram for the 
concatenated data matrix is shown in Fig. 6. As only the ITS 
mutation rate is estimated in previous studies (Leavitt et al. 
2012a, b), asecond chronogram was prepared using this locus 
alone. Results from this analysis have to be treated with cau- 
tion, as the species tree is not strictly clock-like (B. glabra has a 
shorter branch), and the ITS mutation rate has been taken from 
Melanohalea, a lichen-forming genus in the same family. Both 
analyses produced similar values, and the divergence of the 
B. glabra lineage was estimated at 6.9 Mya (95 Yo HPD = 3.5- 
10.8) in the concatenated matrix analysis, and 6.5 Mya (95 Yo 
HPD = 2.2-11.4) in the ITS data alone. The Ko, NA, and WD 
split was estimated at 1.0 Mya (95 Yo HPD - 0.3-2.2) from the 
concatenated matrix and 0.6 Mya (95 % HPD = 0.2—1.5) from 
the ITS data alone. 


Bayesian Skyline Plots (Fig. 5, left) indicate a recent population 
increase in the NA and WD clades. However, the sequences 
contained few informative mutations and the deepest coales- 
cence was reached in around 700000 yr, with no population 
changes detectable further back from this period. Tests of neu- 
trality (Fig. 5, right) are commonly used to support inferences 
from Bayesian Skyline Plots. As indicated by non-significant 
Tajima’s D and Fs results, all sampled groups seem in mutation- 
drift equilibrium, with the exception of the GAPDH locus of the 
NA clade which had a significant negative D value (Fig. 5 bold). 
This could indicate a recent expansion or 'purifying' selection, 
as seen in the concatenated Bayesian Skyline analysis, but 
other loci did not support this hypothesis. 


Markers ITS, IGS, and GAPDH indicate population stability over 
the recent past for clades NAand WD, with putative even more 
recent small population expansions. Due to the low variability 
of the loci, and the putative loss of demographic signals, this 
hypothesis is not confirmed by this analysis. 


Integrated assessment of datasets 


Depending on the analysis, different numbers of putative spe- 
cies were suggested, ranging from four to six (Table 4). All 
analyses, however, confirmed that the combination of morpho- 
logical and chemical characters generally used for species cir- 
cumscription in the complex was inadequate. GAPDH, despite 
its low variability, was the only marker tested that supported 
species-rank assignations for the clades G, Ko, NA, and WD 
(Table 3). ITS, one of the most used loci for DNA barcoding in 
lichen-forming fungi, did not unambiguously distinguish those 
clades. The new markers FRBi15 and FRBi16, despite their 
higher variability, showed inconclusive results and putative 
recombination events. The microsatellite data (Fig. 2) supported 
the DNA sequences results and reflected internal variability not 
revealed in our sequence data, showing that the WD cluster 
was much more diverse than NA, which had a particularly low 
diversity. 
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Fig.5 Bryoria sect. Implexae. — Left. Bayesian Skyline Plots for each clade predicted by the ITS marker and the concatenated loci matrix. The X-axis of each 
graph represent time (in Myr), and the Y-axis represents the value for the log of the effective population size as relative changes, because generation times in 
Bryoria species are unknown. Grey shadows indicate the upper and lower 95 % credible intervals. — Right. Results from neutrality tests for each marker and 
clade, indicating (in bold) any statistically significant deviation from neutrality. — h = number of haplotypes; Fs = Fu’s Fs; D = Tajima’s D. 


Table 4 Summary of the number of putative species suggested by the different methods used for each dataset in Bryoria sect. Implexae. 


Method 


Data 


Figure / reference 


Number of putative species 


Traditional concept DNA sequences and phenotypic Velmala et al. (2014) 12 

Chemical Phenotypic Appendix 2 Left c.4 

Morpho-chemical Phenotypic Appendix 2 Right Not conclusive 

Phylogeny DNA sequences of ITS, IGS, and GAPDH Fig. 6 4=G+Ko+NA+WD 

Phylogeny DNA sequences of FRBi15 Fig. 1 Not conclusive 

Phylogeny DNA sequences of FRBi16 Fig. 1 Not conclusive 

STRUCTURE Microsatellites Fig. 2 5 = G + Ko + NA + WDr + WDg 

PCoA Microsatellites Fig. 3 4 =G + (Ko, NA) + WDr + WDg 
Haplotype Network DNA sequences of ITS, IGS, and GAPDH Appendix 4 4 =G + Ko + NA + WD 

ABGD DNA sequences of ITS, IGS, and GAPDH Table 3 4 =G + Ko + NA + WD 

PTP DNA sequences of ITS, IGS, and GAPDH Table 3 5 = G + Ko + NA + WDr + WDg 
GMYCs DNA sequences of ITS, IGS, and GAPDH Table 3 6 =G + Ko + NA + pik5 + WDr + WDg 
GMYCm DNA sequences of ITS, IGS, and GAPDH Table 3 5 =G + Ko + NA + WDr + WDg 
DISSECT DNA sequences of ITS, IGS, and GAPDH Fig. 4 5= G + Ko + NA + pik5 + WD 
TAXONOMY sent or present, tuberculate or fissural, white to dark. Isidia or 


Bryoria sect. Implexae (Gyeln.) Brodo & D. Hawksw., Opera 
Bot. 42: 114. 1977 


Basionym. Bryopogon sect. Implexae Gyeln., Feddes Repert. Spec. Nov. 
Regni Veg. 38: 223, 238. 1935. 


Type species. Bryoria implexa (Hoffm.) Brodo & D. Hawksw. 1977. = 
Usnea [unranked] implexa Hoffm. 1796. = Bryoria fuscescens (Gyeln.) Brodo 
& D. Hawksw. 1977; but see below. 


Species with a fruticose, hair-like, subpendent to mainly pen- 
dent thallus, lateral spinules or spinulose branches absent, 
whitish grey to brown or black, often paler in the basal parts. 
Angles between branches variable, acute to obtuse or even 
rounded. Pseudocyphellae absent or present, then frequently 
inconspicuous, + fusiform, concolorous or whitish. Soralia ab- 


isidioid spinules absent. Apothecia mainly absent, if present, 
usually afunctional. Chemistry varied, with no detectable or with 
one or a combination of major substances, including alectori- 
alic, barbatolic, connorstictic, fumarprotocetraric, gyrophoric, 
norstictic, protocetraric, psoromic and possibly salazinic acids, 
atranorin, and chloroatranorin. Photobiont Trebouxia ‘hypogym- 
niae’ (Lindgren et al. 2014). 


Notes — Most species included in Brodo & Hawksworth 
(1977) under Bryoria sect. Implexae were transferred to other 
sections in Myllys et al. (2011). In the light of our results (but 
see the Discussion later), Bryoria sect. Implexae includes the 
four species treated below. Comments on particular morpho- 
logical or chemical traits that may be helpful for distinguishing 
these taxa are given under each species. Nevertheless, nearly 
all cited characters are shared by different taxa, so they can 
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be interpreted as ‘cryptic’. The species names adopted here 
are epitypified by sequenced material here in order to fix their 
identities at the molecular level. This epitypification is essential 
to fix the application of these names as no DNA sequences are 
available and cannot be obtained from old type material of most 
species names. The old types cannot therefore be critically 
identified for purposes of the precise application of the names 
and so epitypes may be designated (Turland et al. 2018: Art. 
9.9 and Ex 9). 


As molecular data are necessary for unambiguous species level 
identification in the taxonomy proposed here, we recommend 
using the collective 'Bryoria fuscescens complex' when referring 
to material lacking molecular data. 


Bryoria fuscescens (Gyeln.) Brodo & D. Hawksw., Opera 
Bot. 42: 83. 1977 


Basionym. Alectoria fuscescens Gyeln., Nytt Mag. Naturvidensk. 70: 55. 
1932, nom. cons. (cf. Hawksworth & Jørgensen 2013). 

Synonyms. Lichen chalybeiformis L., Sp. PI. 2: 1155. 1753, (nom. cons.) 
nom. rej. against Bryoria fuscescens (cf. Hawksworth & Jørgensen 2013). 

Bryoria chalybeiformis (L.) Brodo & D. Hawksw., Opera Bot. 42: 81. 1977. 

Usnea [unranked] implexa Hoffm., Deutschl. Fl., Zweiter Teil: 134. 1796. 

Bryoria implexa (Hoffm.) Brodo & D. Hawksw., Opera Bot. 42: 121. 1977. 

Parmelia jubata Q. [var.] capillaris Ach., Methodus, Sectio post.: 273. 
1803. 

Bryoria capillaris (Ach.) Brodo & D. Hawksw., Opera Bot. 42: 115. 1977. 

Alectoria jubata C. [var.] lanestris Ach., Lichenogr. Universalis: 593. 1810. 

Bryoria lanestris (Ach.) Brodo & D. Hawksw., Opera Bot. 70: 88 1977. 

Alectoria kuemmerleana Gyeln., Nytt Mag. Naturvidensk. 70: 49. 1932. 

Bryoria kuemmerleana (Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 
155. 1977. 

Alectoria prolixa var. subcana Nyl. ex Stizenb., Ann. Naturhist. Mus. 
Wien 7: 129. 1892, nom. rej. against Bryoria fuscescens (cf. Hawksworth & 
Jørgensen 2013). 

Bryoria subcana (Nyl. ex Stizenb.) Brodo & D. Hawksw., Opera Bot. 42: 
91. 1977. 

Alectoria vrangiana Gyeln., Magyar Bot. Lapok 31: 46. 1932. 

Bryoria vrangiana (Gyeln.) Brodo & D. Hawksw., Opera Bot. 42: 97. 1977. 


Type specimens. FiNLAND, Tavastia austr., Hollola, ad truncos Pini locis 
apricioribus in silva, Sept. 1882, J.P. Norrlin (Norrlin, Herb. Lich. Fenn. No. 
46) (BP 33947 - lectotype designated by Hawksworth 1972: 217). — FINLAND, 
Etela-Savo, Savitaipale, 600 m NW of Mustapaa, 61, N1721° E27,6900°, 
2005, L. Myllys 464 (HA.H92097 15 (L139)) — epitype designated here, Myco- 
Bank MBT381730. 


Nomenclature — A large number of species rank names 
belong to this group, and are synonymised, but these have not 
been epitypified with sequenced material. Further information 
on synonyms and type materials can be seen in Hawksworth 
(1972), Brodo & Hawksworth (1977) and Velmala et al. (2014). 
Although no samples of Bryoria austromontana have been 
studied, the published description and images (J@rgensen & 
Galloway 1983) suggest this taxon also belong here. 


The earliest species rank epithets amongst these are chalybei- 
formis dating from 1753 (Lichen chalybeiformis), and implexa 
dating from 1799 (Usnea implexa). The former has been re- 
jected against Bryoria fuscescens, but not against other species 
names apart from B. subcana (Hawksworth & Jorgensen 2013). 
A proposal to add the four earlier names Alectoria capillaris, 
Usnea implexa, A. kuemmerleana, and A. lanestris to the two 
against which Alectoria fuscescens is already conserved is 
being prepared separately. Protection against A. vrangiana is 
not required as it appeared in the same work as A. fuscescens. 
While the proposal is under discussion, the name B. fuscescens 
should be adopted in accordance with Rec. 14A.1 (Turland et 
al. 2018). 


We refrained from epitypifying and taking up any of the earlier 
and potentially competing names by epitypification primarily as 
the name B. fuscescens is the most commonly used species 


name in the complex, is well established, the most widely used* 
and is already conserved over two earlier species names in the 
complex. In addition, all the other names have been associated 
with particular morphotypes or chemotypes since the 1970s, 
and so their use might be mistaken as applying to a taxon with 
those particular traits. 


If the proposal for rejection of the previously mentioned compet- 
ing synonyms is not accepted, the principle of priority would 
rule the use of the earliest and not already rejected, validly 
published name at the species rank, i.e., Usnea implexa (and 
then the combination Bryoria implexa), which would require 
epitypification by sequenced material in order to fix the precise 
application of that name. The species was first described from 
Germany but with no named locality, and neotypified by an 
unlocalised and undated specimen in Hoffmann's herbarium 
in Moscow (Hoffmann 8569, MW) which may be part of the 
original material from Germany or have been collected later 
and perhaps in Russia (Hawksworth 1969a). As the neotype 
contained psoromic acid, and the epithet has therefore been 
applied to that chemotype, a potential sequenced epitype 
should represent that chemotype and ideally also have been 
collected in Germany. No such specimen was available to us 
during this study. 

Distribution — Widely distributed, known from cool temper- 
ate to boreal and arctic areas of Europe, Asia, North America, 
and Africa. There are also reports from Antarctica, Oceania, 
and South America, but material from those regions has not 
yet been studied molecularly and so we cannot confirm that 
they belong to this complex. 


Notes — Bryoria fuscescens is highly variable in morphology 
and chemistry, and many of the analysed specimens develop 
soralia. Further, atranorin, which is not normally detectable 
in the other three species accepted here, is frequently found 
in concentrated extracts from both sorediate and esorediate 
morphs. 


Bryoria glabra (Motyka) Brodo & D. Hawksw., Opera Bot. 42: 
86. 1977 


Basionym. Alectoria glabra Motyka, Fragm. Florist. Geobot. 6: 448. 1960. 


Type specimens. USA, Washington, Olympic Peninsula, Clallam Co., Hur- 
ricane Ridge, 5800 ft, on trunk of Abies lasiocarpa, 24 July 1950, B.I. Brown 
& W.C. Muenscher 129 (US — holotype). — USA, Alaska, Mainland, Valley 
between the Bucher and Gilkey Glaciers, southern end of subalpine valley, 
on east side of creek running through valley, subalpine forest, N58°47'20.12" 
W134°30'0.10", 773 m elevation, on Tsuga mertensiana twigs, 4 Aug. 2011, 
K.L. Dillman 4Aug11:1 (UBC (L406)) — epitype designated here, MycoBank 
MBT381731. 


Distribution — Known from northern Europe (Scandinavia), 
and North and South America. In North America, it is most 
abundant in the Pacific North-West. 


Notes — Distinguishing features in well-developed speci- 
mens are the brownish thallus with a regular branching pat- 
tern, generally with obtuse and rounded angles between the 
branches, and broad oval and usually white soralia. It is, how- 
ever, difficult to separate poorly developed or small specimens 
conclusively, so molecular sequences are recommended for 
unambiguous identifications. Only fumarprotocetraric and pro- 
tocetraric acids have been detected in this species, and these 
are characteristically produced in the soralia. 


The Alaskan specimen is selected as the epitype here as se- 
quences are available from all loci, whereas the material we 


* Hits obtained for these names in Google and Google Scholar respectively 
on 12 April 2018 were: B. fuscescens (12300 and 1620), B. capillaris 
(10300 and 999), B. implexa (6020 and 447), B. kuemmerleana (226 and 
20), B. lanestris (6 740 and 185), and B. vrangiana (532 and 61). 
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have seguenced from Washington state (type locality) only has 
data on the ITS locus. 


Bryoria kockiana Velmala, Myllys & Goward in Velmala etal., 
Ann. Bot. Fenn. 51: 361. 2014 


Type specimen. USA, Alaska, Fairbanks, North Star Borough, 26 July 
2011, D. Nossov 20019-1 (UBC (L394) — holotype). 


Distribution — Known only from Alaska (USA) and British 
Columbia (Canada), on conifer branches. 


Notes — Few specimens of this species have so far been 
studied, and these are characterised by the absence of any 
whitish grey tone in the thallus, the lack of soralia, and greyish 
to brown branches with conspicuous, white to concolourous, 
broad, elongate-fusiform, sometimes slightly raised pseudo- 
cyphellae. It lacks TLC-detectable substances or produces 
psoromic acid. The not validly published designation Alectoria 
krogii D. Hawksw. 1972 may be synonymised here. 


Bryoria pseudofuscescens (Gyeln.) Brodo & D. Hawksw., 
Opera Bot. 42: 127. 1977 


Basionym. Alectoria pseudofuscescens Gyeln., Ann. Hist.-Nat. Mus. Natl. 
Hung. 28: 283. 1934, and Rev. Bryol. Lichénol. 7: 51. 1934. 

Synonyms. Bryoria friabilis Brodo & D. Hawksw., Opera Bot 42: 118. 
1977. 

Bryoria pikei Brodo & D. Hawksw., Opera Bot 42: 125. 1977. 

Bryoria inactiva Goward et al., Ann. Bot. Fenn. 51: 360. 2014. 


Type specimens. USA, Oregon, Benton County, Corvallis, on old apple 
trees, Dec. 1931, EP. Sipe 669 (BP 33958 — holotype of Alectoria pseudofusce- 
scens). — Canana, British Columbia, 25 Sept. 2006, T. Goward 07-02-2011 
(UBC (S222) - epitype selected here, MycoBank MBT381732; British Colum- 
bia, Clearwater Valley, 0.5 km S of Philip Creek, 'Edgewood West, 715 m, 
9 Nov. 2011, T. Goward 11-61 (UBC (L347) — holotype of Bryoria inactiva). 


Nomenclature — Anumber of species rank names are syno- 
nymised to this taxon, but these have not been epitypified with 
sequenced material. All these names, however, are later in 
date than pseudofuscescens, and so could not have priority 
over that name. Further information on type materials can be 
seen in Brodo & Hawksworth (1977) and Velmala et al. (2014). 
Although no samples of Bryoria salazinica have been studied at 
molecular level, the published description and images (Brodo 
& Hawksworth 1977) suggests this taxon also belong here. 

Distribution — Only known from North America, growing on 
bark, branches, rock or soil. 


Notes — Characterised by the absence of soralia and detect- 
able atranorin. 


DISCUSSION 


Phylogenetic relationships 


Species concepts in Bryoria sect. Implexae have previously 
been based primarily on well-characterised northern European 
and North American specimens (Hawksworth 1972, Brodo & 
Hawksworth 1977, Velmala et al. 2014). Velmala et al. (2014) 
recognised 11 species on the basis of morphological and 
chemical characters, but many of these were not supported by 
molecular data, and different species names were accepted for 
taxa that could not be distinguished molecularly. We discovered 
that these demarcations broke down when specimens from 
more southern European populations were incorporated. This 
is shown in a phenetic analysis using only phenotypically diag- 
nostic characters (Appendix 2), where the resulting groups are 
not resolved as clear-cut morphospecies. Indeed, any character 
previously used in the group could be used to define the three 
lineages of the Bryoria fuscescens complex (Fig. 6). 


Sexual structures are of major importance in species identifi- 
cation in fungi, but here the rarity of apothecial production 
has hampered their study in most Bryoria species. Any such 
features found would in any case be of limited practical value 
in identification as nearly all samples lack apothecia. Extrolite 
composition has been accorded a major role in species delimi- 
tation in the complex, but many of the substances that were 
considered to be of diagnostic value are biosynthetically closely 
related, being produced by the same gene cluster (pks genes; 
Keller & Hohn 1997), and may be environmentally influenced 
(Myllys et al. 2016, Lutsak et al. 2017). 


Integrative taxonomy, rather than phylogenies based only on 
neutral markers, are increasingly being used to resolve com- 
plex taxonomical groups (e.g., Dayrat 2005, Will et al. 2005, 
Lumley & Sperling 2011, Zamora et al. 2015, Caparrós et al. 
2016). Microsatellites are also now widely used in intraspecific 
population studies because of their high variability (Widmer et 
al. 2012, Dal Grande et al. 2014), and in species complexes 
with diffuse genetic barriers, microsatellite data can improve 
DNA sequence resolution (Lumley & Sperling 2011, Vanhaecke 
et al. 2012). It is generally assumed that DNA sequence data 
reflect the evolution of the species, but these data only reflect 
the history ofthe studied loci, which may sometimes be different 
from the species history overall (Nichols 2001). In this case, 
we demonstrated that traditionally used loci (ITS, IGS, and 
GAPDH) and microsatellite data reveal similar clades, whereas 
other intergenic loci (FRBi15 and FRBi16) produced discrepant 
but statistically supported lineages. These incongruences may 
be due to recombination, hybridisation, or incomplete lineage 
sorting, as documented in many other species groups (eg. 
Jakob & Blattner 2006, McGuire et al. 2007, Edwards et al. 
2008, Stewart et al. 2014). In lichen-forming fungi, outcrossing 
and recombination have been demonstrated, for example, in 
Lobaria pulmonaria (Zoller et al. 1999, Singh et al. 2012, Kel- 
ler & Scheidegger 2016), Letharia (Kroken & Taylor 2001a, b, 
Altermann et al. 2014), and Cladonia (Steinová et al. 2013). 


Apothecia are usually absent in Bryoria sect. Implexae, and 
even when present may not contain mature spores. If cryptic 
sexuality is not occurring, hybridization is unlikely to provide an 
explanation of our data. In the absence of sexual reproduction, 
any recombination is improbable, although some fungi lacking 
sexual structures show recombination events attributable to 
parasexual cycles (Schoustra et al. 2007). We did detect signals 
suggesting putative recombination in the FRBi loci, but not in 
the standard three loci used in the taxonomy adopted here. 
Recombination signals may reflect some mitotic recombination, 
actual or ancient sexual reproduction (Douhan et al. 2007) or be 
merely false positives produced by chance production of similar 
sequences. In any case, recombination alone is insufficient 
to explain all the discordances found. For example, only one 
putative recombination event was detected in FRBi15, and 
disentangling the FRBi16 recombination points is insufficient 
to obtain the topology of the three-locus phylogeny. Incongru- 
ences may also be caused by the analysis of different paralogs 
of FRBi15 and FRBi16 amplified with the new primers, but this 
seems improbable, as no paralogs have been detected in the 
SSRs of these loci (Boluda et al. unpubl.). However, our results 
indicate recent diversification and large effective population 
sizes in this lichenised complex. Thus, incongruences amongst 
loci seem rather attributable to incomplete lineage sorting. 


The different putative species concepts generated by the spe- 
cies delimitation programs (Table 3) may not be only due to 
the different algorithms applied, but also because some of the 
programs were designed for use in single-locus phylogenies 
(e.g., ABGD, PTP, or GMYC). Nevertheless, all the clustering 
analyses showed a tendency to distinguish four groups, G, Ko, 
NA, and WD (Table 4; Fig. 6). STRUCTURE was unable to 
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define these groups until reaching the K = 6 hypothesis, which 
can be attributed to the highly unbalanced sampling sizes; the 
analyses shows that the WD cluster is much more variable 
than NA, which was not divided into subgroups until K = 10. 


Specimen 49 (identified as B. pikei, Fig. 2 asterisk), probably 
reflects the impossibility of unequivocally distinguishing that 
species from B. capillaris; however, as sequence amplification 
of this specimen failed, we cannot determine if this mismatch 
was due to misidentification or DNA contamination. Haplotype 
network analyses have been extensively used in infraspecific 
population and less frequently closely related species groups 
(Houbraken et al. 2012, Pino-Bodas et al. 2016). Even this 
type of nested clade phylogeographic analysis has some critics 
(e.g., Knowles 2008, Templeton 2009). The resulting network 
in the present case is concordant with that obtained from other 
analyses. If two DNA barcoding standard marker networks are 
obtained from a single analysis with a 95 % parsimony con- 
nection limit, members of each network might be considered 
as different species (Hart & Sunday 2007), showing the clear 
isolation of B. glabra from the other taxa in the complex. In the 
case of connections with several mutation steps, as between 
clades Ko, NA, and WD, taxon delimitation is below the species 
level, but in any case indicates some degree of genetic isolation. 


The relationships amongst the Ko, NA, and WD clades remains 
unresolved, indicating that the evolutionary history may be too 
complex to be adequately captured by dichotomous phylo- 
genies based on a few neutral markers. Moreover, the putative 
presence of shared ancestral polymorphisms amongst the 
clades may be producing incompatible topologies, which result 
in clades with low support. 


We also performed analyses to estimate changes in past po- 
pulation sizes, which may have affected current clade diver- 
sity. Genealogies of most plant and animal species coalesce 
between 2.58—0.01 Mya (Grant 2015), and our estimated 
intervals are within this range. However, in our case the dates 
are relatively recent, with the oldest coalescences at 0.7 Mya. 
A flat graphic is generally interpreted as population stability 
but can also be due to a lack of detection power produced by 
small sample sizes. Moreover, a small rise in the curve near the 
present, seen in Fig. 5, can be a consequence of the random 
sampling of the MCMC haplotype trees (Grant 2015); this result 
must therefore be treated with caution. Our sequences do not 
bear imprints of ancient population history, but rather more re- 
cent population growth, for example by extensions northwards 
in post-glacial periods. The loss of information may also arise 
from bottlenecks (i.e., a marked reduction in population size), 
local extinctions, and subsequent recolonization. Additionally, 
the use of genes with low levels of polymorphisms, may impede 
a robust reconstruction of population sizes through time. 


Species concept 


Some species delimitation analyses, such as STRUCTURE, 
GMYC, PTP, can overestimate the number of taxa meriting 
formal recognition, particularly when sampling is uneven or in 
species with a strong intraspecific genetic structure (Altermann 
et al. 2014, Modica et al. 2014, Alors et al. 2016, Del-Prado et 
al. 2016). ABGD, in contrast, has been considered as rather 
conservative, less prone to species overestimation and less 
sensitive to unbalanced sampling. While that method only de- 
tects discontinuities in DNA sequence variation (Puillandre et 
al. 2011), it can also be expected to fail in species with strong 
population genetics structures, for example ones containing 
exclusive haplotypes. All species delimitation programs will 
provide a number of reasonably discrete groups ('evolving line- 
ages’) that should be evaluated for consideration as meriting 
species rank, but the decision has to be by taxonomists with 
experience in the group concerned. Some of our analyses 


suggest that groups such as WDr, WDg, or pik5 might merit 
species rank, but our experience, together with the results from 
other analyses shown here, leads us to reject this hypothesis. 
We conclude that the most pragmatic solution, supported by 
the general trend of the results from the different analyses we 
performed is to consider clades G, Ko, NA, and WD as the 
species Bryoria glabra, B. kockiana, B. pseudofuscescens, 
and B. fuscescens, respectively. 


Clade age can contribute to decisions as to species limits. 
Divergence time estimates can be robust if the analyses are 
performed with well-resolved phylogenies and can incorporate 
fossil calibration points, as in some vertebrates (Perelman et 
al. 2011). Contrarily, in lichenised fungi, fossils are rare and in 
many cases enigmatic or with ambiguous relationships (Thomas 
et al. 2014, Hawksworth 2015, Kaasalainen et al. 2015). In ad- 
dition, generation times can be expected to be different among 
species, as is the case with nulTS locus mutation rates between 
herbaceous and woody plants, or even a difference of almost 
an order of magnitude between different plant genera (Kay 
et al. 2006). Here we used a nulTS mutation rate estimated 
from Melanohalea, a genus in the same family (Leavitt et al. 
2012a), species of which frequently grow with Bryoria and 
reproduce asexually as well as sexually. The split of B. glabra 
from the other taxa in the B. fuscescens complex is estimated 
at c. 6.9 Mya, and clearly separated from the much later 
divergence of the other three species estimated at c. 1 Mya 
(0.6 Mya if only the nulTS locus is used). This contrasts with 
other lichenised species considered of recent origin, estimated 
around 2.5—5 Mya (Pliocene), with any Pleistocene speciation 
event rare and always older than 1 Mya (Amo de Paz et al. 
2012, Leavitt et al. 2012a, b, Molina et al. 2016). As the three 
B. fuscescens complex clades seem to have diverged more 
recently, the extent of their reproductive isolation is uncertain, 
and the discovery of intermediate lineages amongst other 
named species from unsampled geographical regions, such 
as continental Asia remains possible. 


In the absence of supporting phenotypic, geographic, or eco- 
logical differences, the recent divergence, and the possibility 
of incomplete lineage sorting, clades Ko, NA, and WD may be 
considered as conspecific evolving lineages. It is, however, 
important to recognise the lineages formally in order to facilitate 
their conservation by enabling their threat status to be assessed 
by IUCN criteria. We decided not to adopt the rank of subspe- 
cies as that is now hardly used in mycology, and then not in 
any consistent way; traditionally this rank was used in plants for 
morphologically distinguishable populations with geographical 
differences and where intermediates occurred where they were 
sympatric (Stuessy 2009). 

The formal recognition of cryptic lineages at species level, as 
suggested by our analyses, emerges as the most appropriate 
solution. Cryptic speciation is now recognised as a common 
phenomenon in Parmeliaceae, and our results are in accord 
with other studies in which molecular markers in combination 
with statistical tools revealed genetically distinct lineages pre- 
viously hidden under a single taxon name in this family (e.g., 
Singh et al. 2015, Alors et al. 2016, Del-Prado et al. 2016, 
Divakar et al. 2016, Leavitt et al. 2016). Further, this solution 
is in line with the increasing need to formally recognize cryp- 
tic species-level lineages in all fungi (Hibbett 2016); indeed, 
cryptic speciation may mean that there are on average ten or 
more fungi masked in formally named species (Hawksworth & 
Lücking 2017). 

Of the lineages recognised here, only the WD clade emerged 
as cosmopolitan, occurring in Europe, Asia, North America, and 
Africa (Appendix 5). NA and Ko have been collected so far only 
in North America, despite our extensive sampling in Europe 
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(Appendix 5). Further sampling, especially in South America, 
Asia, and Africa, is needed before any finer-scale biogeographic 
patterns might be detected. 


The practical issue of naming older museum specimens and 
material in ecological surveys could be resolved by recognis- 
ing the three groups as species within a broad concept, such 
as an aggregate, complex, or adding ‘s.lat.’. We considered 
commending the adoption of the suffix ‘agg.’ for material when 
precise molecular species identifications cannot be made. While 
this has been done in a few other groups of fungi (e.g., Parn- 
men et al. 2013, Pazoutova et al. 2015), ‘complex’ has come to 
be used more widely and was strongly favoured at the Cryptic 
Speciation in Fungi symposium in Utrecht in September 2017 
(report awaited). We therefore suggest the use of ‘complex’ here 
but recognise some may prefer to use ‘agg.’ or 's.lat.'. Where 
DNA samples can be obtained and analysed, we recommend 
use of the GAPDH locus, as all the other tested markers are 
not able to distinguish with confidence the three species we 
recognize in the complex. 


Infraspecific phenotypic variation 


While our results support rejection of the morphospecies con- 

cept in this group of lichens, two main phenotypes can nev- 

ertheless often be distinguished by the naked eye in the field: 

1. the pale grey ‘capillaris’ morphotype (including B. capillaris 
and B. pikei, Fig. 6b); and 

2. the fuscous brown to dark brown 'fuscescens' morphotype, 

in which most other species names are placed (Fig. 6a, c). 


The chemical characters are not always checked by field 
workers, and while the 'capillaris' morphotype typically has 
benzyldepsides, the 'fuscescens' morphotype lacks those 
compounds and has fumarprotocetraric acid or various depsi- 
dones. However, there are dark morphs with benzyldepsides 
(once called f. fuscidula), and pale grey ones with fumarpro- 
tocetraric acid (e.g., B. subcana) or other depsidones (e.g., B. 
kuemmerleana). It is conceivable that the two morphotypes 
originated before the separation of B. pseudofuscescens and 
B. fuscescens, as both colour variants and chemistries appear 
in both taxa. This phenomenon cannot be explained by a simple 
ongoing speciation event in which one lineage has originated 
new adaptations, but is still not isolated from the parental 
lineages, as neither are monophyletic in a paraphyletic clade. 


The difference in phenotype cannot be attributed to different 
algal partners as all material in the complex shares the same 
species and even in many cases the same nulTS haplotypes of 
Trebouxia (Lindgren et al. 2014, Boluda et al. unpubl.). Further, 
as we used neutral markers to detect gene-flow gaps between 
lineages, the phenotypes are also not the result of genetic isola- 
tion, and other possibilities must be considered. 


Ithas recently been reported that yeast morphs ofthe lichenicol- 
ous and gall-forming basidiomycete genus Cyphobasidium can 
be abundant in or on the outermost cortical tissues of Bryoria 
species (Spribille et al. 2016). Spribille et al. (2016) reported 
a possible relation between Cyphobasidium yeast abundance 
and vulpinic acid production in two other species of Bryoria, 
B. fremontii and B. tortuosa, and also visualised these yeasts 
in material identified as B. capillaris phenotypes. Contrary to 
the claims of Spribille et al. (2016), these fungi do not appear 
to be an integral part of the mutualism (Oberwinkler 2017). It 
is, however, feasible that the yeasts cells are able to develop 
to a greater extent in ‘capillaris’ morphotypes as the cortices 
can have lumpier polysaccharide deposits than do those of 
‘fuscescens’ (Hawksworth 1969b, Boluda et al. 2014, Esseen 
et al. 2017). How the occurrence of yeast morphs of these 
lichenicolous fungi in the surface of the cortex could possibly 
determine colour morphotypes is obscure. 


Material referred to the 'capillaris and 'fuscescens' pheno- 
types has been reported to show slight differences in water 
holding capacity (Esseen et al. 2015), and also the pigments 
may provide protection against excesses of light (Fárber et al. 
2014). Further, in southern Europe particularly, the 'capillaris' 
phenotype tends to be restricted to humid, shaded, and pro- 
tected or undisturbed environments than the 'fuscescens' one, 
something already recognised by Motyka (1964). Additionally, 
in northern Europe, dark specimens containing barbatolic acid 
are much more common than in southern Europe, where they 
are extremely rare (Myllys et al. 2016, Esseen et al. 2017). As 
both phenotypes can grow side by side and even intermixed 
on the same trees, where environmental conditions must be 
identical, ecological plasticity has to be discounted. However, 
some unknown epigenetic modification could perhaps have a 
role in that process, as once a metabolic pathway is activated 
or silenced, it may be hardly modifiable under more or less 
neutral environmental conditions, transferring the phenotypes 
to the clonal offspring. Specimens with dark thalli, containing 
barbatolic acid, or with pale thalli with traces of barbatolic and 
also containing other extrolites, could represent transitional 
specimens. 


Molecular and morphological rates of divergence may some- 
times be uncoupled. Incomplete lineage sorting arises when 
anancestral polymorphism persists through a speciation event 
and each polymorphism can lead to different alleles being car- 
ried among descendants (Maddison 1997, Hartl & Clark 2007). 
Consequently, different tree topologies may be obtained de- 
pending which specimens or loci are used. Rosenberg (2003) 
has shown that 5.3N, generations are needed for a species to 
acquire monophyly at 99 Yo of its loci given that all loci in the 
sister species are also monophyletic. That indicates that for a 
species of 1 M individuals with a generation time of 10 yr, the 
full monophyly will only be reached 50 M years after speciation, 
whereas only around 1000 yr may be needed for species with 
small populations (Naciri & Linder 2015). Incomplete lineage 
sorting may be frequent in closely related taxa or during a spe- 
ciation process (Hobolth et al. 2011, Blanco-Pastor et al. 2012, 
Saag et al. 2014, Naciri & Linder 2015), as may be considered 
in our case. The topological incongruence observed among the 
standard loci, FRBi15 and FRBi16, supports the incomplete line- 
age sorting hypothesis as one of the main reasons explaining 
why morphospecies are not monophyletic. While neutral mark- 
ers are useful for understanding gene-flow patterns, adaptive 
markers provide the evolutionary pressure that contributes to 
speciation (Emelianov et al. 2004, Hey 2006, Holderegger et al. 
2006). As adaptive markers are under natural selection, certain 
alleles can be present in some morphospecies and absent in 
others, even if there is gene flow amongst them (Lumley & 
Sperling 2011). The use of phylogenomic datasets may provide 
a more accurate and supported phylogenetic reconstruction, 
especially if the appropriate scale of loci variability is selected 
from all the genome (Leavitt et al. 2016). However, if there 
are high levels of incomplete lineage sorting, it might not be 
expected that morphospecies would appear forming supported 
clades. Nevertheless, genomic data may reveal few mutations 
linked to certain morphospecies, which would be producing 
adaptive traits. Darwin's finches are an iconic example of a rapid 
speciation process, in which there is a mismatch between the 
phylogenetic species concept and phenotype-based taxonomy; 
in that case, genomic studies have detected specific loci sub- 
jected to selection pressure, which are directly related with the 
development of taxon-specific phenotypes (Lamichhaney et al. 
2015). In Bryoria, supposed adaptive traits may be influenced 
by the genes involved in the production of certain extrolites or 
in the epicortical substances (Boluda et al. 2014), which may 
produce differential selection pressure for each morphotype, 
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at least in some environments. The process might be similar 
to that of natural selection of the pale and melanic morphs 
of the Peppered Moth (Biston betularia) in Europe (Majerus 
2009), impeding the fixation of a single morphotype in all 
populations. In our case also, high levels of incomplete lineage 
sorting mixed with a few phenotypically important genes under 
variable degrees of selection in different environments, may 
explain the mismatch observed between phenotypes and geno- 


types. 
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Appendix 1 Chemical, geographical and morphological characters used of Bryoria sect. Implexae samples in the phenogram reconstruction (Fig. 1). 
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Appendix 2 Phenograms based on a presence/absence distance matrix in Bryoria sect. Implexae from: a. Extrolite composition alone; b. extrolite composition, 
with geographical, and morphological data. — Bold branches represent supported clades (bootstrap = 70 Yo, approximately unbiased 2 95 Yo). — Ale. = Alec- 
torialic acid; Bar. = Barbatolic acid; Fum. = Fumarprotocetraric acid; Gyr. = Gyrophoric acid; No subs. = No substances detected; Nor. = Norstictic acid; 


Pso. = Psoromic acid; * = Except specimen named Bryoria capillaris L14.02. 
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Appendix 3 Microsatellite fragment lengths of Bryoria sect. Implexae analysed specimens. 


Sample Bi01 Bi03 Bi04 Bi05 Bi10 Bi11 Bi12 Bi14 Bi19 
capillaris_L01-17 103 279 323 128 437 316 100 365 346 
capillaris __L06-10 103 279 323 128 437 316 100 365 346 
capillaris __L07-15 103 279 323 128 437 316 100 365 346 
capillaris L08-12 103 279 323 128 434 316 103 361 346 
capillaris_L13-03 94 279 325 138 434 318 115 361 346 
capillaris_L14-02 112 279 327 128 437 320 100 365 346 
capillaris L141 123 281 323 136 434 316 100 361 346 
capillaris L15-15 112 279 323 - 437 316 100 365 346 
capillaris_L16-21 103 279 323 128 437 316 100 365 346 
capillaris_L211 112 279 323 136 434 316 103 361 346 
capillaris_L270 112 279 323 138 434 316 124 361 346 
capillaris_S192 112 279 323 128 437 316 100 365 346 
capillaris_S2 94 279 323 138 434 316 124 361 346 
friabilis 02 114 281 317 137 434 310 100 - 350 
friabilis L355 120 277 316 137 436 310 100 365 350 
friabilis L407 132 281 316 137 434 310 131 365 346 
friabilis S395 109 277 316 132 438 310 100 365 350 
fuscescens_L12-03 112 279 323 138 434 316 103 361 352 
fuscescens_L12-05 94 281 323 138 434 316 103 361 350 
fuscescens_L139 94 279 323 138 434 316 103 361 352 
fuscescens_L149 123 279 323 138 434 316 103 361 350 
fuscescens_L15-21 123 277 325 138 434 318 115 361 352 
fuscescens_L160 94 281 323 136 434 316 103 361 352 
fuscescens_L189 112 279 323 138 434 316 100 361 352 
fuscescens_L224 112 279 323 138 434 316 103 361 352 
fuscescens_L232 94 279 323 138 434 316 103 361 352 
fuscescens_L305 117 279 323 138 434 316 137 361 352 
fuscescens_S109 112 281 323 138 434 316 118 361 350 
fuscescens_S157 117 279 323 136 434 316 103 361 350 
fuscescens_S24 112 281 323 138 434 316 106 361 352 
fuscescens_S256 94 281 323 - 434 316 124 363 354 
fuscescens_S259 94 277 323 138 437 316 103 363 352 
fuscescens_S260a 82 279 323 138 437 316 103 363 352 
fuscescens_S261 82 279 323 138 437 316 103 363 352 
fuscescens_S267 82 279 323 138 437 316 103 363 352 
fuscescens_S272 82 279 323 138 437 316 103 363 352 
fuscescens_S274 94 279 323 138 434 316 - 361 350 
fuscescens_S369 - 279 323 138 437 316 103 363 352 
fuscescens_S379 94 279 323 138 437 316 100 363 352 
fuscescens_S380 112 279 323 138 434 316 118 361 352 
fuscescens_S56 123 279 323 136 434 316 100 361 350 
glabra 01 114 283 306 132 426 299 115 369 344 
glabra 02 120 283 306 132 426 299 - 369 344 
glabra_03 120 283 306 132 426 299 - 369 344 
glabra_04 120 283 306 132 426 299 - 369 344 
glabra 05 120 283 306 132 426 299 - 369 344 
glabra_L186 114 283 304 132 426 297 115 371 344 
glabra_L406 114 283 304 132 426 297 115 371 344 
glabra_L414 114 283 306 132 426 299 115 371 344 
glabra_S388 114 283 306 132 426 299 115 371 344 
implexa_L01-01 135 281 323 138 434 316 106 361 350 
implexa L06-05 112 263 323 136 434 316 103 361 350 
implexa L10-03 106 279 323 128 437 316 100 365 346 
implexa L11-15 94 279 323 136 435 316 103 361 352 
implexa L16-15 112 281 325 138 434 318 115 361 352 
implexa S168 112 279 323 138 434 316 115 361 350 
implexa_S22 94 279 323 136 436 316 118 361 352 
implexa_S36 94 279 325 136 434 318 103 361 346 
implexa_S39 - - 323 138 434 316 - 361 350 
implexa_S67 94 279 325 136 434 318 103 361 346 
inactiva_L206 114 277 317 137 434 311 100 365 350 
inactiva_L323b 114 281 316 132 434 310 131 365 350 
inactiva L347 114 277 317 132 434 310 124 365 350 
inactiva L358 109 281 317 132 436 310 106 365 350 
inactiva_S239a 109 277 314 132 434 308 100 365 350 
inactiva_S384 114 277 316 137 434 310 100 365 350 
inactiva_S392 120 281 316 137 434 310 100 365 350 
kockiana_L394 94 279 317 136 472 310 109 365 344 
kockiana_L396 94 279 317 136 472 310 109 365 344 
kuemmerleana L04-03 129 279 325 138 434 318 103 361 352 
kuemmerleana L09-04 112 281 323 128 437 316 100 365 346 
kuemmerleana L09-07 112 281 323 128 437 316 100 365 346 
kuemmerleana L16-17 112 281 323 138 434 316 106 361 352 
kuemmerleana L244 117 279 323 138 434 316 115 361 350 
kuemmerleana L274 94 279 323 136 434 316 103 361 352 
kuemmerleana_L275 94 279 323 136 434 316 103 361 352 
kuemmerleana_S128 100 279 323 136 434 316 100 361 354 


kuemmerleana_S160 117 279 323 138 434 316 103 361 350 
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Appendix 3 (cont.) 


Sample Bi01 Bi03 Bi04 Bi05 Bi10 Bi11 Bi12 Bi14 Bi19 
pikei 02 117 277 317 137 434 310 100 365 344 
pikei 04 117 277 317 137 434 310 100 365 344 
pikei 05 117 277 317 137 434 310 100 365 344 
pikei 07 117 281 317 137 434 311 109 - 344 
pikei 09 114 281 317 132 - 310 100 - 344 
pikei_10 117 277 317 137 434 311 100 365 344 
pikei_11 114 277 317 137 434 310 100 365 344 
pikei 12 114 277 - 137 434 - 100 - 344 
pikei_13 - 277 317 137 434 310 118 365 344 
pikei_14 117 277 317 137 434 310 127 365 344 
pikei 15 114 281 317 137 436 310 109 - 344 
pikei_a 114 - 323 132 437 316 100 - 346 
pikei_b 117 277 317 137 434 311 100 365 344 
pikei_c 108 - - 137 - 314 100 326 344 
pikei_d 117 277 - 137 437 310 100 - 344 
pikei_L197 114 277 316 132 434 310 109 365 344 
pikei_L210 109 277 316 132 434 310 100 365 344 
pikei_L241 126 277 316 137 434 310 109 - 352 
pikei_L374 114 277 316 137 434 310 106 365 344 
pikei_L376 132 277 316 137 434 310 127 365 352 
pikei_L377 109 277 316 137 434 310 103 365 344 
pikei_S221 109 277 316 132 436 310 106 365 352 
pikei_S362 114 281 317 137 434 311 100 365 344 
pikei_S368 112 281 316 132 436 310 109 365 344 
pikei_S382 114 281 317 137 434 311 100 365 344 
pikei_S383a 114 281 317 132 436 311 100 365 344 
pikei_S390 114 277 316 137 434 310 100 365 344 
pikei_S394 126 277 316 137 434 310 100 365 352 
pseudofuscescens_S222 - 277 316 132 - 310 100 365 - 
pseudofuscescens_S232 114 277 317 137 436 310 100 365 350 
pseudofuscescens_S370 117 277 316 132 434 310 100 365 350 
pseudofuscescens_S371 117 281 316 132 434 310 100 365 350 
pseudofuscescens_S377 - 277 - 132 434 311 100 365 - 
pseudofuscescens_S386 112 277 317 132 438 311 100 365 350 
pseudofuscescens_S387 117 277 316 132 434 310 127 365 350 
sp L395 94 273 317 136 460 310 109 365 344 
sp S392 97 283 317 136 472 310 118 365 344 
vrangiana L02-20 112 279 323 138 434 316 106 361 352 
vrangiana_L03-07 94 279 325 136 434 318 115 365 346 
vrangiana L05-17 117 281 323 138 434 316 121 361 350 
vrangiana L07-03 94 281 323 138 434 316 118 361 352 
vrangiana L07-19 123 279 323 138 434 316 103 361 350 
vrangiana_L08-19 94 281 323 136 436 316 103 361 352 
vrangiana_L08-20 94 279 323 138 434 316 121 361 352 
vrangiana L10-13 123 283 323 128 437 316 100 365 346 
vrangiana L12-11 94 279 323 138 434 316 144 361 352 
vrangiana L13-12 94 279 323 138 434 316 103 361 352 
vrangiana_L272 112 279 323 136 434 316 - 361 350 
vrangiana_L273 94 279 323 138 436 316 103 361 350 
vrangiana_L300 94 281 323 136 434 316 115 361 350 
vrangiana_L307 123 279 325 138 434 318 103 361 352 
vrangiana_S10 123 279 323 136 434 316 124 361 352 
vrangiana_S164 94 279 323 136 434 316 127 361 350 
vrangiana_S166 117 279 323 138 434 316 144 365 352 
vrangiana_S196a 94 281 323 138 434 316 118 361 350 
vrangiana S341 - 279 323 - 434 316 - 361 350 
vrangiana_S385 94 279 323 138 434 316 103 361 350 
vrangiana_S42 117 279 323 138 434 316 131 361 352 
vrangiana_S45 112 279 323 136 434 316 100 361 350 
vrangiana_S57 94 279 323 138 434 316 115 361 352 
vrangiana_S59 94 279 323 138 434 316 115 361 352 
vrangiana_S6 123 279 323 138 434 316 103 361 350 
vrangiana_S62 112 279 323 138 434 316 115 361 352 
vrangiana_S72 94 279 323 136 436 316 103 361 352 
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Appendix 4 Haplotype network in Bryoria sect. Implexae of a concatenated matrix containing ITS, IGS and GAPDH sequences. The analysis coded gaps 
as missing data and used a 95 Yo connection limit. Numbers represent the specimens shown in Table 1 and colours depict the STRUCTURE microsatellites 
genepool (Fig. 2). Connecting line length do not depict the genetic distance. Each line represents a single mutation connected by black small circles. Circle 
size is related with the number of analysed specimens. — * = WDb (Wide Distributed brown cluster) specimens. 
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Appendix 5 Distribution of Bryoria sect. Implexae specimens examined. Two samples of geographical interest, not analysed in this study, have been added: 
Bryoria kockiana (Canada, British Columbia, 1982, Goward 82-1040, UBC — paratype; cf. Velmala et al. 2014) and Bryoria fuscescens (Tanzania, Kilimanjaro, 
2016, Boluda & Kitara, MAF-Lich. 20750). — Basemap source: U.S. National Park Service (NPS) Natural Earth physical map. 


